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We describe a newly developed cosmological hydrodynamics code based on the 
weighted essentially non-oscillatory (WENO) schemes for hyperbolic conservation 
laws. High order finite difference WENO schemes are designed for problems with 
piecewise smooth solutions containing discontinuities, and have been successful in 
applications for problems involving both shocks and complicated smooth solution 
structures. We couple hydrodynamics based on the WENO scheme with standard 
Poisson solver - particle-mesh (PM) algorithm for evolving the self-gravitating 
system. A third order total variation diminishing (TVD) Runge-Kutta scheme 
has been used for time-integration of the system. We brief the implementation 
of numerical technique. The cosmological applications in simulating intcrgalactic 
medium and Lyo forest in the CDM scenario are also presented. 

1 Introduction 

Due to the highly non-Hnearity of gravitational clustering in the universe, one signif- 
icant feature emerging in cosmological hydrodynamic flow is extremely supersonic 
motion around the density peaks developed by gravitational instability, which leads 
to strong shock discontinuities within complex smooth structures. It poses more 
challenge than the typical hydrodynamic simulation without self-gravity. To ad- 
dress this issue, two algorithms have been implemented for high resolution shock 
capturing in cosmology: total-variation diminishing (TVD) scheme^^l and piecewise 
parabohc method (PPM)Pl. 

An alternative hydrodynamic solver to discretize the convection terms in the 
Euler equations is the fifth order finite difference WENO (weighted essentially non- 
oscillatory) method'^l, with a low storage nonlinearly stable third order Runge- 
Kutta time discretization!^]. The WENO schemes, first constructed in third order 
finite volume version by Liu et al.'^l, are based on the essentially non-oscillatory 
scheme (ENO) developed by Harten et al.I^] in the form of finite volume scheme 
for hyperbolic conservative laws. The ENO scheme generalizes the total variation 
diminishing (TVD) scheme of Harten. The TVD schemes typically degenerate to 
first-order accuracy at locations with smooth extrema while the ENO and WENO 
schemes maintain high order accuracy there even in multi-dimensions. The WENO 
method is very robust and stable for solutions containing strong shocks and complex 
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solution structures'^! '[^1. It uses the idea of adaptive stencils in the reconstruction 
procedure based on the local smoothness of the numerical solution to automatically 
achieve high order accuracy and non-oscillatory property near discontinuities. This 
is achieved by using a convex combination of a few candidate stencils, each being 
assigned a nonlinear weight which depends on the local smoothness of the numerical 
solution based on that stencil. WENO schemes can simultaneously provide a high 
order resolution for the smooth part of the solution, and a sharp, monotone shock 
or contact discontinuity transition. 

In the context of cosmological applications, we have developed a hybrid N- 
body/hydrodynamical code that incorporates a Lagrangian particle-mesh algorithm 
to evolve the collision-less matter with the fifth order WENO scheme to solve the 
equation of gas dynamics. This paper describes this code briefly. 

2 Numerical Techniques 

2.1 The Basic Equations 

The hydrodynamic equations for baryons in the expanding universe, without any 
viscous and thermal conductivity terms, can be written as, 

U + di¥\U]= f{t,U) (1) 

where the abbreviation di = d/dx^ has been used, x* denote the proper coordinates, 
t/ is a five-component column vector, U = (p, pvi, pv2, pv^, E), p is the comoving 
density, v = {vi} is the proper peculiar velocity, E is the total energy per unit 
comoving volume including both kinetic and internal energy, P is comoving pressure, 
related to the total energy E hy E = P/{j — 1) + |/5v^ where we assume an ideal 
gas equation of state, P = (7 — l)e, where e is the internal energy density and 7 is 
the ratio of the specific heats of the baryon. The left hand side of Eq.(l) is written 
in the conservative form for mass, momentum and energy, the "force" source term 
on right hand side includes the contributions from the expansion of the universe, 
the gravitation. A„et in the energy equation represents the net energy loss due to 
the radiative heating-cooling of the baryonic gas. 

fit, U) = (0, ~pv + pg, -2^E + pv • g - Acooi) (2) 

where g is the peculiar acceleration in the gravitational field produced by both the 

dark matter and the baryonic matter, which is obtained by solving the Poisson 
equation using standard particle-mesh technique. 

2.2 Hydrodynamic Solver: Finite Difference WENO Schemes 

The fifth order WENO finite diff'erence spatial discretization to a conservation law 
such as 

Ut + f{u)x + g{u)y = (3) 

approximates the derivatives, for example f{u)x, by a conservative difference 
f{u)x\x=xj « (/?+i/2 - /j-1/2) 
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where fj+i/i is the numerical flux. g{u)y is approximated in the same way. Hence 
finite difference methods have the same format for one and several space dimensions, 
which is a major advantage. For the simplest case; of a scalar equation (3) and if 
— the fifth order finite difference WENO scheme has the flux given by 

where three third order accurate fluxes on three different stencils given 

by 

/i+1/2 + ^/K) + 

/]+\/2 = \f{Uj) + - \f{Uj+2)- 

Notice that the combined stencil for the flux ./j+1/2 is biased to the left, which is 
upwinding for the positive wind direction due to the assumption f'{u) > 0. The 
key ingredient for the success of WENO scheme relies on the design of the nonlinear 
weights Wi, which are given by 

Wi _ 7/c 

where the linear weights 7/5 are chosen to yield fifth order accuracy when combining 
three third order accurate fluxes, and are given by 

_ 1 _ 3 _ 3 

71 - 72 - 73 - 

the smoothness indicators arc given by 

/3i = (/(t.,_2) - 2/(u,_i) + /(«,))^ + i (/(«,_2) - 4/(t.,_i) + Zf(ui)f 
/32 = § {liu^-x) - 2f{uj) + f{uj+,)f + \ ifiuj.^) - f{uj+r)f 

Ps = ^ ifiuj) - 2/(«,+i) + .fiu,+2) f + I (3/K) - 4/(«,+i) + .f{u,+2) f , 

and they measure how smooth the approximation based on a specific stencil is in 
the target cell. Finally, e is a parameter to avoid the denominator to become zero 
and is usually taken as £ = 10~^ in the computation. 

This finishes the description of the fifth order finite difference WENO schemel^l 
in the simplest case. As we can see, the algorithm is actually quite simple and the 
user does not need to tunc any parameters in the scheme. 

The WENO finite difference scheme is characterized by the following properties. 
(1) The scheme is proven to be uniformly fifth order accurate including at smooth 
extrema, and this is verified numerically. (2) Near discontinuities the scheme pro- 
duces sharp and non-oscillatory discontinuity transition. (3) The approximation 
is self-similar. That is, when fully discretized with the Runge-Kutta methods dis- 
cussed in next section §2.3, the scheme is invariant when the spatial and time 
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variables are scaled by the same factor. This is a major advantage for approxi- 
mating conservation laws which are invariant under such scaling. For the details of 
how the scheme can be generalized in a more complex situation, eventually to 3D 
systems such as Euler equations, we refer to the reference [7]. 



2.3 Time Discretizations 

To further discretize in time, we use a class of high order nonlinearly stable Runge- 

Kutta time discretizations. A distinctive feature of this class of time discretizations 
is that they are convex combinations of first order forward Euler steps, hence they 
maintain strong stability properties in any semi-norm (total variation norm, maxi- 
mum norm, entropy condition, etc.) of the forward Euler step, with a time step re- 
striction proportional to that for the forward Euler step to be stable, this proportion 
coefficient being termed CFL coefficient of the high order Runge-Kutta method'^l . 
The most popular scheme in this class is the following third order Runge-Kutta 
method for solving ut = L{u,t), where L{u,t) is a spatial discretization operator: 

= Ati(u",t") 

= Jw" + + i AfL(w(i) , r + At) (4) 

o o o ^ 

which is nonlinearly stable with a CFL coefficient 1.0. However, for our purpose 

of 3D calculations, storage is a paramount consideration. We thus use the third 
order low storage nonlinearly stable Runge-Kutta method which was proven to be 
nonlinearly stable in [4], with a CFL coefficient 0.32. 

The timestep is chosen by the minimum value among three time scales. The 
first is from the Courant condition given by 

A.< , ,^^^,xf^- , ^ (5) 
max(|ux| + Cs, \Uy\ + c^, |Wz| -|- Cg) 

where Ax is the cell size, Cg is the local sound speed, Ux, Uy and Uz are the local 
fluid velocities and CFL is the Courant number, typically, we take CFL ~ 0.3. The 
second constraint is imposed by cosmic expansion which requires that Aa/a < 0.02 
within single time step. The last constraint comes from the requirement that a 
particle move no more than a fixed fraction of cell size in one time step. 



3 Cosmological Application: IGM and Lya Forest at High Redshift 

We run our hybrid N-body/hydrodynamic code to compute the cosmic evolu- 
tion of coupled system of both dark matter and baryonic matter in a flat low 
density CDM model (ACDM), which is specifled by the cosmological parame- 
ters (fim, f^A, /i, CTg, rJb) = (0.3,0.7,0.65,0.9,0.035). The simulations were per- 
formed in a periodic, cubical box of size 12h~^Mpc with a 128"^ grid and an 
equal number of dark matter particles. Atomic processes including ionization, 
radiative cooling and heating are modeled as in Cenl^l in a primeval plasma of 
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Figure 1. Scatter plot of IGM temperature vs. density relation randomly drawn from a ACDM 
simulation at redshift z = 3.0. The line is the least-square-fitting of power-law using each set of 
points displayed in the figure. 



hydrogen and helium of composition {X — 0.76, Y — 0.24). The uniform UV- 
background of ionizing photons is assumed to have a power-law spectrum of the 
form J(y) = J21 x 10~^^(i^/i^/f/)~"ergs~"'^cm~^sr~^Hz~-'^ parameterized by pho- 
toionizing flux J21 at the Lyman limit frequency vm^ and is suddenly switched on 
at z 6 to heat the gas and reionize the universe. 

The evolution of the IGM temperature relies on the reionization history of hydro- 
gen and helium, and plays important role in our understanding the Lya absorption 
spectrum. It is expected that a power-law relation between IGM temperature and 
density T = Tq{p/ < p >)'''~-'^ could be estabhshed due to the ionizing photon heat- 
ing the gas'^^l. We test T versus p relation using the simulation samples. Fig.(l) 
shows the scatter plot of temperature-density relation drawn from the simulation 
at the output of redshift z=3. Apparently, the T-p relation approximates roughly 
a power law with a best-fitting value of 7 = 0.56. The dispersions in temperature 
at a fixed density mainly arise from the shock-heating. 

Mock Lya spectra are produced along random lines of sight in the simulation 
box, and meanwhile, the corresponding one-dimension distributions of IGM density 
and peculiar velocity are also extracted from the sample. The results are demon- 
strated in Fig. (2). To compare with the observation, e.g. the Keck spectrum of 
QSO q0014-|-8118, we Gaussian smooth the spectra to match with the spectral res- 
olutions of observation and normalize the spectra by the observed mean flux decre- 
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Figure 2. An example of simulated Lya absorption spectrum and corresponding IGM physical 
properties along randomly selected line of sight. Upper panel: transmitted flux; middle panel: 
density distribution of the neutral hydrogen; lower panel: peculiar velocity distribution. 

ment. For q0014+8118, we have D = 1- < e~'^ >w 0.31. It could be shown that 
there are good agreements of flux power spectrum and some higher-order statistical 
behaviors'^^1 (e.g. intermittency) between the observation and simulations'^^]. 
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